/*
 * ConvexHull.h
 *
 * Created 6/21/2009 By Johnny Huynh
 *
 * Version 00.00.01 6/21/2009
 *
 * Copyright Information:
 * All content copyright © 2008 Johnny Huynh. All rights reserved.
 */

 #ifndef CONVEXHULL_H
 #define CONVEXHULL_H
 
 #include "global.h"
 
 template <typename TYPENAME> class ConvexHull;
 
 #include <vector>
 #include <list>
 #include "OpenGL_Headers.h"
 #include "GL_Triangle.h"
 #include "introsort.h"
 
 typedef ConvexHull<GLfloat> ConvexHullf;
 typedef ConvexHull<GLdouble> ConvexHulld;
 
 /**
  * Class prototype for ConvexHull
  */
 template <typename TYPENAME>
 class ConvexHull
 {
    // Data Members
    //protected:   
    public: 
       std::vector< GL_Triangle<TYPENAME> > triangles; // the set of triangles that form this convex hull
    
    // Local Functions
    public:
        ConvexHull( const GL_Triangle<TYPENAME>* const triangles_Ptr, size_t size );
        ConvexHull( const GL_Triangle<TYPENAME>* const * const triangle_Ptrs_Ptr, size_t size );
        ConvexHull( const ConvexHull<TYPENAME>& ch );
        ~ConvexHull( void );
        inline void operator=( const ConvexHull<TYPENAME>& ch );
    
    // Private Functions
    private:
        inline void quick_hull( Vector3<TYPENAME>* uniqueVertices, const size_t size );
        
    // Private Friend Functions
    // Should probably move these 3 functions to be public in Vector3.h
        template <typename TYPENAME>
        friend inline bool compareByXAxis( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b );
        template <typename TYPENAME>
        friend inline bool compareByYAxis( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b );
        template <typename TYPENAME>
        friend inline bool compareByZAxis( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b );
        template <typename TYPENAME>
        friend inline bool compareByYXZ( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b );
    
    // Friend Functions
    public:
        template <typename TYPENAME> friend inline void render( const ConvexHull<TYPENAME>& ch );
    
    // Private Classes
        template <typename TYPENAME>
        class VertexAboveFacet
        {
        // Data Members
        public:
            const Vector3<TYPENAME>* vertex_Ptr; // pointer to the vertex
            TYPENAME projection; // the magnitude of the vertex from the facet projected onto the normal of the facet
            
        // Local Functions
            VertexAboveFacet( const Vector3<TYPENAME>* const vertex_Ptr = NULL, const TYPENAME projection = 0.0f ) 
                              : vertex_Ptr( vertex_Ptr ), projection( projection ) {}
            VertexAboveFacet( const VertexAboveFacet<TYPENAME>& vaf ) 
                              : vertex_Ptr( vaf.vertex_Ptr ), projection( vaf.projection ) {}
            ~VertexAboveFacet() {}
            inline VertexAboveFacet<TYPENAME>& operator=( const VertexAboveFacet<TYPENAME>& vaf ) 
            { 
                vertex_Ptr = vaf.vertex_Ptr;
                projection = vaf.projection;
                return *this;
            }
            
        // Friend Functions
            template <typename TYPENAME>
            friend inline bool operator<( const VertexAboveFacet<TYPENAME>& va, const VertexAboveFacet<TYPENAME>& vb )
                                                                             { return va.projection > vb.projection; }
        };
        template <typename TYPENAME>
        class Facet // i.e. triangle in 3D space
        {
        // Data Members
        public:
            GL_Triangle<TYPENAME> triangle; // triangle representing this facet
            std::vector< const VertexAboveFacet<TYPENAME> > verticesAbove; // list of vertices above this facet
            std::list< Facet<TYPENAME>* > neighbor_Ptrs; // list of facet neighbors (as pointers) of this facet
                                                                 // two facets are neighbors if they share a ridge (i.e. edge)
        
        // Local Functions    
            Facet( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b, const Vector3<TYPENAME>& c ) :triangle( a, b, c ){}
            Facet( const GL_Triangle<TYPENAME>& triangle ) : triangle( triangle ) {}
            Facet( const Facet<TYPENAME>& facet ) : triangle( facet.triangle ), 
                                                verticesAbove( facet.verticesAbove ), neighbor_Ptrs( facet.neighbor_Ptrs ) {}
            ~Facet() {}
            // A vertex, v, is above/outside a facet if for any vertex on the facet, fv, the line v-fv, projected onto
            // the normal vector, u, of the facet is positive (i.e. dotProduct( (v - fv), u ) > 0 ).
            inline TYPENAME getVertexAboveProjection( const Vector3<TYPENAME>& v ) const
                                    { return (v - triangle.getEndPointA()) * normal( triangle ); }
        // Friend Functions
            template <typename TYPENAME>
            friend bool inline areNeighbors( const Facet<TYPENAME>& f, const Facet<TYPENAME>& g )
            {   // Two facets are neighbors if they share an edge
                Vector3<TYPENAME> fa( f.triangle.getEndPointA() );
                Vector3<TYPENAME> fb( f.triangle.getEndPointB() );
                Vector3<TYPENAME> fc( f.triangle.getEndPointC() );
                
                Vector3<TYPENAME> ga( g.triangle.getEndPointA() );
                Vector3<TYPENAME> gb( g.triangle.getEndPointB() );
                Vector3<TYPENAME> gc( g.triangle.getEndPointC() );
                
                if ( fa == ga )
                {
                    return fb == gb || fb == gc || fc == gb || fc == gc;
                }
                else if ( fa == gb )
                {
                    return fb == ga || fb == gc || fc == ga || fc == gc;
                }
                else if ( fa == gc )
                {
                    return fb == gb || fb == ga || fc == gb || fc == ga;
                }
                
                if ( fb == ga )
                {
                    return fa == gb || fa == gc || fc == gb || fc == gc;
                }
                else if ( fb == gb )
                {
                    return fa == ga || fa == gc || fc == ga || fc == gc;
                }
                else if ( fb == gc )
                {
                    return fa == gb || fa == ga || fc == gb || fc == ga;
                }
                
                /* // Unnecessary
                if ( fc == ga )
                {
                    return fa == gb || fa == gc || fb == gb || fb == gc;
                }
                else if ( fc == gb )
                {
                    return fa == ga || fa == gc || fb == ga || fb == gc;
                }
                else if ( fc == gc )
                {
                    return fa == gb || fa == ga || fb == gb || fb == ga;
                }*/
                
                return false;
            }
        };
        template <typename TYPENAME>
        class Ridge // edge in 3D space
        {
        // Data members
        public:
            const Vector3<TYPENAME>& a; // an end point of the ridge
            const Vector3<TYPENAME>& b; // another end point of the ridge
        
        // Local Functions
            Ridge( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b ) : a(a), b(b) {}
            Ridge( const Ridge<TYPENAME>& r ) : a(r.a), b(r.b) {}
            ~Ridge() {}
            inline Ridge<TYPENAME>& operator=( const Ridge<TYPENAME>& r )
            {
                memcpy( this, &r, sizeof( Ridge<TYPENAME> ) );
                return *this;
            }
        // Friend Functions
            template <typename TYPENAME> friend inline bool operator==( const Ridge<TYPENAME>& r, const Ridge<TYPENAME>& s )
            {
                return ( (r.b == s.a && r.a == s.b) || (r.a == s.a && r.b == s.b) );
            }
        };
 };
 
 /** LOCAL FUNCTIONS **/
 
 /**
  * Constructor - takes in the array of triangles and number triangles that constitute the object
  * we are computing the convex hull around.
  */
 template <typename TYPENAME>
 ConvexHull<TYPENAME>::ConvexHull( const GL_Triangle<TYPENAME>* const triangles_Ptr, size_t size )
 {
    // the object does not exist
    if ( size == 0 )
        return;
    
    // L = set of (unique) unprocessed points of the object for finding the convex hull
    
    // set of unique points of the object
    std::vector< Vector3<TYPENAME> > uniqueVertices;
    getUniqueVertices( uniqueVertices, triangles_Ptr, size );
    
    // the object is a point
    if ( uniqueVertices.size() == 1 )
    {
        triangles.push_back( triangles_Ptr[0] );
        return;
    }
    
    // Perform the quick hull algorithm
    quick_hull( &uniqueVertices[0], uniqueVertices.size() );
 }
 
 /**
  * Alternative Constructor
  */
 template <typename TYPENAME>
 ConvexHull<TYPENAME>::ConvexHull( const GL_Triangle<TYPENAME>* const * const triangle_Ptrs_Ptr, size_t size )
 {
    // the object does not exist
    if ( size == 0 )
        return;
    
    // L = set of (unique) unprocessed points of the object for finding the convex hull
    
    // set of unique points of the object
    std::vector< Vector3<TYPENAME> > uniqueVertices;
    getUniqueVertices( uniqueVertices, triangle_Ptrs_Ptr, size );
    
    // the object is a point
    if ( uniqueVertices.size() == 1 )
    {
        triangles.push_back( *triangle_Ptrs_Ptr[0] );
        return;
    }
    
    // Perform the quick hull algorithm
    quick_hull( &uniqueVertices[0], uniqueVertices.size() );
 }
 
 /**
  * Alternative Constructor
  */
 template <typename TYPENAME>
 ConvexHull<TYPENAME>::ConvexHull( const ConvexHull<TYPENAME>& ch )
                    : triangles( ch.triangles )
 {
    
 }

 /**
  * Destructor
  */
 template <typename TYPENAME>
 ConvexHull<TYPENAME>::~ConvexHull( void )
 {
    
 }
 
 /**
  * operator=() copies the content of the specified ConvexHull
  * to this ConvexHull.
  *
  * @param (const ConvexHull<TYPENAME>&) ch
  */
 template <typename TYPENAME>
 inline void ConvexHull<TYPENAME>::operator=( const ConvexHull<TYPENAME>& ch )
 {
    if ( ch.triangles.size() > 0 )
        triangles.assign( &ch.triangles[0], &ch.triangles[ ch.triangles.size() - 1 ] + 1 );
 }
 
 /** PRIVATE FUNCTIONS **/
 
 /**
  * quick_hull() performs the quick hull algorithm to calculate the convex hull for the specified set of unique vertices.
  * The triangles constructing the convex hull are stored in the set of triangles for this ConvexHull class.
  *
  * @param (const Vector3<TYPENAME>*) uniqueVertices - unique vertices pointer
  * @param (const size_t) size - number of unique vertices
  */
 template <typename TYPENAME>
 inline void ConvexHull<TYPENAME>::quick_hull( Vector3<TYPENAME>* uniqueVertices, const size_t size )
 {
    size_t numberOfUniqueVertices( size );
    
    // the object is a line
    if ( numberOfUniqueVertices == 2 )
    {
        triangles.push_back( GL_Triangle<TYPENAME>( uniqueVertices[0], uniqueVertices[1], uniqueVertices[0] ) );
        return;
    }
    // the object is a triangle
    else if ( numberOfUniqueVertices == 3 )
    {
        triangles.push_back( GL_Triangle<TYPENAME>( uniqueVertices[0], uniqueVertices[1], uniqueVertices[2] ) );
        return;
    }
    
    // 1. Find top 4 points (i.e. ordered on the y-axis) from L and remove the points from L.
    
    // sort the vertices by the Y-Axis placing vertices with a greater y-axis value at the front.
    introsort< Vector3<TYPENAME> >( uniqueVertices, numberOfUniqueVertices, compareByYAxis );
    
    // 2. Construct a simplex (i.e. tetrahedron) from the 4 points as follows:
	//          a) Select 4 points and call them A, B, C, D
	//          b) Pretend we have a triangle A, B, C
	//          c) Find the normal vector of triangle A, B, C as follows: (B - A) x (C - A)
	//          d) Call dotProduct( (D - A), normal_vector ) the D_Direction, 
	//             which is the relative direction of D to the normal vector of the triangle
	//          e) if ( D_Direction > 0 ), the simplex is constructed of the following facets:
	//              	- A -> B -> D	// FA
	//              	- A -> C -> B	// FB
	//              	- A -> D -> C	// FC
	//              	- B -> C -> D	// FD
	//             else // D_Direction <= 0, and the simples is constructed of the following facets:
	//              	- A -> D -> B	// FA
	//              	- A -> B -> C	// FB
	//              	- A -> C -> D	// FC
	//              	- B -> D -> C	// FD
    
    // 3. Add each facet to a list, M.
    
    
    // list M
    std::list< Facet<TYPENAME> > facets;
    
    // Build the Tetrahedron
    {
        Vector3<TYPENAME>& a( uniqueVertices[0] );
        Vector3<TYPENAME>& b( uniqueVertices[1] );
        Vector3<TYPENAME>& c( uniqueVertices[2] );
        Vector3<TYPENAME>& d( uniqueVertices[3] );
    
        // Pretend we have a triangle a, b, c
        Vector3<TYPENAME> normal_vector( crossProduct( (b - a), (c - a) ) ); // normal vector of triangle a, b, c
        TYPENAME d_direction( (d - a) * normal_vector ); // the direction of d relative to the normal of triangle a, b, c
        
        if ( d_direction > 0 )
        {
            // A -> B -> D	// FA
            facets.push_back( Facet<TYPENAME>( a, b, d ) );
	        // A -> C -> B	// FB
	        facets.push_back( Facet<TYPENAME>( a, c, b ) );
	        // A -> D -> C	// FC
	        facets.push_back( Facet<TYPENAME>( a, d, c ) );
	        // B -> C -> D	// FD
	        facets.push_back( Facet<TYPENAME>( b, c, d ) );
	    }
	    else // ( d_direction <= 0 )
	    {
	        // A -> D -> B	// FA
            facets.push_back( Facet<TYPENAME>( a, d, b ) );
	        // A -> B -> C	// FB
	        facets.push_back( Facet<TYPENAME>( a, b, c ) );
	        // A -> C -> D	// FC
	        facets.push_back( Facet<TYPENAME>( a, c, d ) );
	        // B -> D -> C	// FD
	        facets.push_back( Facet<TYPENAME>( b, d, c ) );
	    }
    }
    
    size_t uniqueVertex_index = 3;
    
    //    A point is above/outside a facet if for any point on the facet, fp, the line p-fp, projected onto
    //    the normal vector, u, of the facet is positive (i.e. dotProduct( (p - fp), u ) > 0 ).
    // 4. For each point in L, if the point is above/outside a facet, link the point to the facet. Each facet should
    //    contain a list of point pointers, which indicates the points that are above the facet.
    // for each unique vertex
    
    // for each unique vertex
    {
        size_t vertex_index( numberOfUniqueVertices );
        if ( vertex_index > uniqueVertex_index )
            do
            {
                --vertex_index;
                
                // for each facet
                std::list< Facet<TYPENAME> >::iterator facet_itr( facets.begin() );
                if ( facet_itr != facets.end() )
                    do
                    {
                        // link the vertex to the facet if the vertex is above the facet
                        TYPENAME m( (*facet_itr).getVertexAboveProjection( uniqueVertices[ vertex_index ] ) );
                        if ( m > EPSILON )
                        {
                            (*facet_itr).verticesAbove.push_back( VertexAboveFacet<TYPENAME>( 
                                                                                &uniqueVertices[ vertex_index ], m ) );
                        }
                        ++facet_itr;
                    } while ( facet_itr != facets.end() );
                
            } while ( vertex_index > uniqueVertex_index );
    }
    
    // 6. For each facet, f, in M,
	//    if facet, f, does not have an outside point,
	//    remove f from M and insert f into N.
    {
        // for each facet
        std::list< Facet<TYPENAME> >::iterator facet_itr( facets.end() );
        if ( facet_itr != facets.begin() )
            do
            {
                --facet_itr;
                
                // if facet, f, does not have an above/outside point,
	            // remove f from M and insert f into N.
                if ( (*facet_itr).verticesAbove.empty() )
                {
                    triangles.push_back( (*facet_itr).triangle );
                    // delete the facets that have no above vertices
                    facet_itr = facets.erase( facet_itr );
                }
            } while ( facet_itr != facets.begin() );
    }
    
    // A facet is a neighbor of a specific facet, f, if the facet shares an edge with f.
    // 7. Link the facets to its neighbors that has at least one outside point; each facet should contain a list 
    //    of pointers to its neighbor facets that have at least one outside point. Do this as follows:
    //    For each facet, f, in M,
	//      add all other facets in M to be f's neighbors list
	{
        // for each facet
        std::list< Facet<TYPENAME> >::iterator facet_itr( facets.begin() );
        if ( facet_itr != facets.end() ) // != facets.end() - 1 )
            do
            {
                // for each other facet
                std::list< Facet<TYPENAME> >::iterator facet_itr2( facet_itr );
                ++facet_itr2;
                if ( facet_itr2 != facets.end() )
                    do
                    {
                        (*facet_itr).neighbor_Ptrs.push_back( &*facet_itr2 );
	                    (*facet_itr2).neighbor_Ptrs.push_back( &*facet_itr );
	                    
                        ++facet_itr2;
                    } while ( facet_itr2 != facets.end() );
                    
                ++facet_itr;
            } while ( facet_itr != facets.end() ); //!= facets.end() - 1 );
	}
	
	//int loop_counter( ch_counter );
	
	// 8. While M is not empty (such that M contain at least one facet),
	while ( facets.size() != 0x0 )
	{
	    /*if ( --loop_counter < 0 )
	    {
	        if ( showM )
	        {
	            triangles.clear();
	            while ( !facets.empty() )
	            {
	                triangles.push_back( facets.front().triangle );
	                facets.pop_front();
	            }
	        }
	        break;
	    }*/
	//      find the point, p, that is furthest in the direction of normal vector of f 
	//      (i.e. max( (pi - fp), u ), and remove it from f.
	    Facet<TYPENAME>& f = facets.front(); //facets[0];
	    introsort< VertexAboveFacet<TYPENAME> >( &f.verticesAbove[0], f.verticesAbove.size() );
	    const Vector3<TYPENAME>& p( *f.verticesAbove[0].vertex_Ptr );
	//  	Create a facet pointer list, V, and initialize V with f.
	    // V is a list of facet pointers where the point, p, is above/outside them
	    std::vector< Facet<TYPENAME>* > V;
	    V.push_back( &f );
	//  	Create an empty facet pointer list, W.   
	    std::vector< Facet<TYPENAME>* > W;
	    
	    {
	//      Create a facet pointer list, neighbor_Ptrs, and initialize it with the neighbor facets of f.
	        std::vector< Facet<TYPENAME>* > neighbor_Ptrs( f.neighbor_Ptrs.begin(), f.neighbor_Ptrs.end() );
	        
	//      For each facet, n, in neighbor_Ptrs
	        // we use a for-loop here because neighbor_Ptrs will dynamically grow in size
	        for ( size_t i = 0; i < neighbor_Ptrs.size(); ++i )
	        {
	            Facet<TYPENAME>& n( *neighbor_Ptrs[ i ] );
	            
	//  		if the point, p, lies outside of n,            
	            if ( n.getVertexAboveProjection( p ) > EPSILON )
	            {
	//              add the neighbors of n to neighbor_Ptrs if the neighbors are not already in neighbor_Ptrs or V
	                // for each neighbor of n
	                std::list< Facet<TYPENAME>* >::iterator n_neighbor_Ptr_itr( n.neighbor_Ptrs.begin() );
	                
	                while ( n_neighbor_Ptr_itr != n.neighbor_Ptrs.end() ) 
	                {
	                    bool isUnvisitedNeighbor( true );
                        // for each neighbor in neighbor_Ptrs
                        size_t neighbor_index( neighbor_Ptrs.size() );
                        // if ( neighbor_index > 0 )    // obviously true
                        do
                        {
                            --neighbor_index;
                            
                            // if the neighbor of n is a visited neighbor
                            if ( *n_neighbor_Ptr_itr == neighbor_Ptrs[ neighbor_index ] )
                            {
                                isUnvisitedNeighbor = false;
                                break;
                            }
                            
                        } while ( neighbor_index > 0 );
                        
                        if ( isUnvisitedNeighbor )
                        {
                            // for each facet in V
                            size_t v_index( V.size() );
                            // if ( v_index > 0 ) // obviously true
                            do
                            {
                                --v_index;
                                
                                // if the neighbor of n is a visited neighbor
                                if ( *n_neighbor_Ptr_itr == V[ v_index ] )
                                {
                                    isUnvisitedNeighbor = false;
                                    break;
                                }
	                        
                            } while ( v_index > 0 );
                            
                            if ( isUnvisitedNeighbor )
                            {
                                // add the neighbor of n to neighbor_Ptrs
                                neighbor_Ptrs.push_back( *n_neighbor_Ptr_itr );
                            }
                        }
                        
                        ++n_neighbor_Ptr_itr;
	                }
	                   
	//  			add n to V    
	                V.push_back( &n );
                }
	//  		else
                else
                {
	//  			add n to W      
                    W.push_back( &n ); 
                }           
            }
	    }
	
	//  	For each facet, v, in V,
	    {
	        size_t v_facet_index( V.size() );
	        // if ( v_facet_index != 0 ) // obviously true
            do
            {
                --v_facet_index;
	//  		For each facet, w, in W,
                size_t w_facet_index( W.size() );
                if ( w_facet_index == 0 )
                    break;
	            else // ( w_facet_index != 0 )
                    do
                    {
                        --w_facet_index;
                        
	                    // for each neighbor of w
                        Facet<TYPENAME>& w( *W[ w_facet_index ] );
                        std::list< Facet<TYPENAME>* >::iterator neighbor_Ptr_itr( w.neighbor_Ptrs.end() );
                        while ( neighbor_Ptr_itr != w.neighbor_Ptrs.begin() )
                        {
                            --neighbor_Ptr_itr;
                            
                            // if v is a neighbor of w
                            if ( V[ v_facet_index ] == *neighbor_Ptr_itr )
                            {
	//  			remove v as a neighbor from w.
                                neighbor_Ptr_itr = w.neighbor_Ptrs.erase( neighbor_Ptr_itr );
                            } 
                        }
                    } while ( w_facet_index > 0 );
            } while ( v_facet_index > 0 );
	    }
	    
	//  	Create an empty list of ridges (i.e. edges), R, which will be the set of horizon ridges.
	    std::vector< Ridge<TYPENAME> > R;
	//  	For each facet in V,
	    {
	        std::vector< Ridge<TYPENAME> > ridges; // set of unique ridges
	        std::vector< Ridge<TYPENAME> > duplicate_ridges;
	        // for each facet in V
	        size_t facet_index( V.size() );
	        if ( facet_index > 0 )
	            do
	            {
	                --facet_index;
	//    		// For a triangle with end points a, b, c, the ridges are a->b, b->c, and c->a.
	                Facet<TYPENAME>& f( *V[ facet_index ] );
	                // We can assume that the 3 ridges are unique compared to one another because each facet is
	                // constructed from unique points (from above).
	                Ridge<TYPENAME> arrayOfRidges[3] = { 
	                                        Ridge<TYPENAME>( f.triangle.getEndPointA(), f.triangle.getEndPointB() ),
	                                        Ridge<TYPENAME>( f.triangle.getEndPointB(), f.triangle.getEndPointC() ),
	                                        Ridge<TYPENAME>( f.triangle.getEndPointC(), f.triangle.getEndPointA() ) };
	                
	                bool isUnique0( true );
	                bool isUnique1( true );
	                bool isUnique2( true );
	//  		For each ridge, e, of the facet in V,	                
	                size_t ridges_index( ridges.size() );
	                if ( ridges_index > 0 )
	                    do
	                    {
	                        --ridges_index;
	                        if ( arrayOfRidges[0] == ridges[ridges_index] )
	                            isUnique0 = false;
	                        if ( arrayOfRidges[1] == ridges[ridges_index] )
	                            isUnique1 = false;
	                        if ( arrayOfRidges[2] == ridges[ridges_index] )
	                            isUnique2 = false;
	                    } while ( ridges_index > 0 );
	                 
	//  			add e to R if e is not shared by any other facet in V (i.e. e is a horizon ridge in V and outlines the 
    //              shape constructed by the set of facets in V)
    //              Notice if two triangles share an edge, the edge must be connected at the same two end points as we 
    //              constructed the basis of the convex hull that way. Also, only up to two triangles can share an edge.
	                if ( isUnique0 )
	                    ridges.push_back( arrayOfRidges[0] );
	                else
	                    duplicate_ridges.push_back( arrayOfRidges[0] );
	                if ( isUnique1 )
	                    ridges.push_back( arrayOfRidges[1] );
	                else
	                    duplicate_ridges.push_back( arrayOfRidges[1] );
	                if ( isUnique2 )
	                    ridges.push_back( arrayOfRidges[2] );
	                else
	                    duplicate_ridges.push_back( arrayOfRidges[2] );
	                
	            } while ( facet_index > 0 );
	            
	        // for each ridge in ridges
	        size_t ridges_index( ridges.size() );
	        if ( ridges_index > 0 )
	            do
	            {
	                --ridges_index;
	                
	                // for each duplicate ridge in ridges
	                size_t duplicate_ridges_index( duplicate_ridges.size() );
	                bool isHorizonRidge( true );
	                if ( duplicate_ridges_index > 0 )
	                    do
	                    {
	                        --duplicate_ridges_index;
	                        
	                        if ( ridges[ ridges_index ] == duplicate_ridges[ duplicate_ridges_index ] )
	                        {
	                            isHorizonRidge = false;
	                            break;
	                        }
	                    } while ( duplicate_ridges_index > 0 );
	                // if the ridge in ridges is a horizon ridge, add the ridge to R
	                if ( isHorizonRidge )
	                    R.push_back( ridges[ ridges_index ] );
	                
	            } while ( ridges_index > 0 );
	    }
	//      Create a numerical variable, numberOfFacetsWithAboveVertices, to keep track of the number of new facets added to M
	    size_t numberOfFacetsWithAboveVertices( 0 );
	    
	
	//  	For each horizon ridge, e, in R    
	    {
	        size_t ridge_index = R.size();
	        if ( ridge_index > 0 )
	            do
	            {
	                --ridge_index;
	//  		create a new facet using the ridge and the point, p (e.g. a->b->p).                
	                Facet<TYPENAME> newFacet( R[ ridge_index ].a, R[ ridge_index ].b, p );
	                
	//  	    for each facet in V
	                size_t v_index( V.size() );
	                if ( v_index > 0 )
	                    do
	                    {
	                        --v_index;
	                        
	//  			for each remaining outside point of the facet in V,
	                        std::vector< const VertexAboveFacet<TYPENAME> >& verticesAbove( V[ v_index ]->verticesAbove );
	                        size_t vertex_index( verticesAbove.size() );                   
	                        if ( vertex_index > 0 )
	                            do
	                            {
	                                --vertex_index;
	                                
                                    VertexAboveFacet<TYPENAME> vertex( verticesAbove[ vertex_index ].vertex_Ptr, 
                                             newFacet.getVertexAboveProjection( *verticesAbove[ vertex_index ].vertex_Ptr ) );
 	//  				if the point is above/outside the new facet, link the point to the new facet
                                    if ( vertex.projection > EPSILON )
                                    {
                                        newFacet.verticesAbove.push_back( vertex );
                                    }
                                    
	                            } while ( vertex_index > 0 ); 
	                    } while ( v_index > 0 );
	                
	//  		if the new facet does not contain at least one outside point
	//  			add the facet to N
	                if ( newFacet.verticesAbove.empty() )
	                    triangles.push_back( newFacet.triangle );
	//  		else // the new facet contain at least one outside point
	                else
	                {
	//          add the new facet to M and refer to this new facet in M by reference for the algorithm below
	                    facets.push_back( newFacet );
	                    Facet<TYPENAME>& newFacetRef( facets.back() );//facets[ facets.size() - 1 ] );
    //          increment numberOfFacetsWithAboveVertices
                        ++numberOfFacetsWithAboveVertices;
	//          For each facet, w, in W,
	                    size_t w_index( W.size() );
	                    if ( w_index > 0 )
	                        do
	                        {
	                            --w_index;
	                            
	                            Facet<TYPENAME>& w( *W[ w_index ] );
	//  				    if w is a neighbor of the new facet
                                Ridge<TYPENAME> newFacetRidge( newFacetRef.triangle.getEndPointA(), 
                                                                newFacetRef.triangle.getEndPointB() );
                                
                                Ridge<TYPENAME> ridgeA( w.triangle.getEndPointA(),
                                                            w.triangle.getEndPointB() );
                                Ridge<TYPENAME> ridgeB( w.triangle.getEndPointB(),
                                                            w.triangle.getEndPointC() );
                                Ridge<TYPENAME> ridgeC( w.triangle.getEndPointC(),
                                                            w.triangle.getEndPointA() );
                                
                                if ( newFacetRidge == ridgeA || newFacetRidge == ridgeB 
                                                             || newFacetRidge == ridgeC )
                                {
    //  					    link w and the new facet as a neighbors
                                    newFacetRef.neighbor_Ptrs.push_back( &w );
                                    w.neighbor_Ptrs.push_back( &newFacetRef );
                                }

                            } while ( w_index > 0 );   
	/*//  		    for each neighbor facet of the facets in V that is not a facet in V,
	                    // for each facet in V
	                    size_t v_index( V.size() );
	                    // if ( v_index > 0 ) // obviously true
                        do
                        {
                            --v_index;
                            Facet<TYPENAME>& facet( *V[ v_index ] );
                            // for each neighbor of the facet
                            size_t neighbor_index( facet.neighbor_Ptrs.size() );
                            if ( neighbor_index > 0 )
                                do
                                {
                                    --neighbor_index;
                                    Facet<TYPENAME>& neighbor( *facet.neighbor_Ptrs[ neighbor_index ] );
                                    
                                    // for each facet in V
                                    bool isNotInV( true );
                                    size_t v_index2( V.size() );
                                    // if ( v_index2 > 0 ) // obviously true
                                    do
                                    {
                                        --v_index2;
                                        
                                        // if the neighbor is a facet in V
                                        if ( &neighbor == V[ v_index2 ] )
                                        {
                                            isNotInV = false;
                                            break;
                                        }
                                        
                                    } while ( v_index2 > 0 );
                                        
                                    // if the neighbor is not a facet in V
                                    if ( isNotInV )
                                    {
    //  				    if the neighbor facet is a neighbor of the new facet // or automatically link them as neighbors for co-planar issue fix. Then again, this is equivalent to saying all facets in M are neighbors of one another.
                                        Ridge<TYPENAME> newFacetRidge( newFacetRef.triangle.getEndPointA(), 
                                                                        newFacetRef.triangle.getEndPointB() );
                                        
                                        Ridge<TYPENAME> ridgeA( neighbor.triangle.getEndPointA(),
                                                                    neighbor.triangle.getEndPointB() );
                                        Ridge<TYPENAME> ridgeB( neighbor.triangle.getEndPointB(),
                                                                    neighbor.triangle.getEndPointC() );
                                        Ridge<TYPENAME> ridgeC( neighbor.triangle.getEndPointC(),
                                                                    neighbor.triangle.getEndPointA() );
                                        
                                        if ( newFacetRidge == ridgeA || newFacetRidge == ridgeB 
                                                                     || newFacetRidge == ridgeC )
                                        {
    //  					    link the facet and the new facet as a neighbors
                                            newFacetRef.neighbor_Ptrs.push_back( &neighbor );
                                            neighbor.neighbor_Ptrs.push_back( &newFacetRef );
    //  					    add new facet as a pointer to newFacet_Ptrs
                                            newFacet_Ptrs.push_back( &newFacetRef );
                                        }
                                        //if ( areNeighbors( newFacet, neighbor ) )
                                        //{
    //  					    link the facet and the new facet as a neighbors
                                            //newFacetRef.neighbor_Ptrs.push_back( &neighbor );
                                            //neighbor.neighbor_Ptrs.push_back( &newFacetRef );
    //  					    add new facet as a pointer to newFacet_Ptrs
                                            //newFacet_Ptrs.push_back( &newFacetRef );
                                        //}
                                    }
                                    
                                } while ( neighbor_index > 0 );
                            
                        } while ( v_index > 0 );*/

	                }
	                    
	            } while ( ridge_index > 0 );
	    }
    //      Use numberOfFacetsWithAboveVertices to access the new facets in M. 
    //      Link these new facets as neighbors to one another if they share an edge.
        {
            // for each new facets in M
            std::list< Facet<TYPENAME> >::iterator facet_itr( facets.end() );
	        while ( numberOfFacetsWithAboveVertices > 0 )
	        {
                --facet_itr;
                --numberOfFacetsWithAboveVertices;
	            
                // for each of the other new facets in M
                size_t numberOfFacetsWithAboveVertices2( numberOfFacetsWithAboveVertices );
                std::list< Facet<TYPENAME> >::iterator facet_itr2( facet_itr );
                while ( numberOfFacetsWithAboveVertices2 > 0 )
                {
                    --facet_itr2;
                    --numberOfFacetsWithAboveVertices2;
                    // if the two new facets are neighbors, link them as neighbors
                    if ( areNeighbors( *facet_itr, *facet_itr2 ) )
                    {
                        (*facet_itr).neighbor_Ptrs.push_back( &*facet_itr2 );
                        (*facet_itr2).neighbor_Ptrs.push_back( &*facet_itr );
                    }
                 }
	         }
	     }
	/*//  	For each facet in newFacet_Ptrs
	    size_t facet_index( newFacet_Ptrs.size() );
	    if ( facet_index > 1 )
	        do
	        {
	            --facet_index;
	            
	//  		for each of the other facets in newFacet_Ptrs
	            size_t facet_index2( facet_index );
	            if ( facet_index2 > 0 )
	                do
	                {
	                    --facet_index2;
	                    
	//  			if the two facets are neighbors, link them as neighbors
	                    if ( areNeighbors( *newFacet_Ptrs[ facet_index ], *newFacet_Ptrs[ facet_index2 ] ) )
	                    {
	                        newFacet_Ptrs[ facet_index ]->neighbor_Ptrs.push_back( newFacet_Ptrs[ facet_index2 ] );
	                        newFacet_Ptrs[ facet_index2 ]->neighbor_Ptrs.push_back( newFacet_Ptrs[ facet_index ] );
	                    }
	                    
	                } while ( facet_index2 > 0 );
    	        
	        } while ( facet_index > 1 );*/
	        
	//  	For each facet in V,
	//  		remove the facet from M
        {
            size_t v_index( V.size() );
            
            if ( v_index > 0 )
                do
                {
                    --v_index;
                    
                    // for each facet in M
                    std::list< Facet<TYPENAME> >::iterator facet_itr( facets.end() );
                    if ( facet_itr != facets.begin() )
                        do
                        {
                            --facet_itr;
	                        
                            // if the facet in M is in V
                            if ( &*facet_itr == V[ v_index ] )
                            {
                                // remove the facet from M
                                facet_itr = facets.erase( facet_itr );
                                break;
                            }
	                    
                        } while ( facet_itr != facets.begin() );
                    
                } while ( v_index > 0 );
        }
        
	}
 }
 
 /** PRIVATE FRIEND FUNCTIONS **/
 
 /**
  * compareByXAxis() returns true if the x-axis value of the specified a is greater than the
  * x-axis value of the specified b; otherwise, false is returned.
  *
  * @param (const Vector3<TYPENAME>&) a
  * @param (const Vector3<TYPENAME>&) b
  */
 template <typename TYPENAME>
 inline bool compareByXAxis( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b )
 {
    return a.x > b.x;
 }
 
 /**
  * compareByYAxis() returns true if the y-axis value of the specified a is greater than the
  * y-axis value of the specified b; otherwise, false is returned.
  *
  * @param (const Vector3<TYPENAME>&) a
  * @param (const Vector3<TYPENAME>&) b
  */
 template <typename TYPENAME>
 inline bool compareByYAxis( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b )
 {
    return a.y > b.y;
 }
 
 /**
  * compareByZAxis() returns true if the z-axis value of the specified a is greater than the
  * z-axis value of the specified b; otherwise, false is returned.
  *
  * @param (const Vector3<TYPENAME>&) a
  * @param (const Vector3<TYPENAME>&) b
  */
 template <typename TYPENAME>
 inline bool compareByZAxis( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b )
 {
    return a.z > b.z;
 }
 
 /**
  * compareByYXZ() returns true if the y-axis value of the specified a is greater than the
  * y-axis value of the specified b. If the y-axis values are equal, then the x-axis values are compared. 
  * Finally, if the x-axis values are equal, then the z-axis values are compared; otherwise, false is returned.
  *
  * @param (const Vector3<TYPENAME>&) a
  * @param (const Vector3<TYPENAME>&) b
  */
 template <typename TYPENAME>
 inline bool compareByYXZ( const Vector3<TYPENAME>& a, const Vector3<TYPENAME>& b )
 {
    if ( a.y == b.y )
    {
        if ( a.x == b.x )
        {
            return a.z > b.z;
        }
        else
            return a.x > b.x;
    }
    else
        return a.y > b.y;
 }
 
 /** FRIEND FUNCTIONS **/
 
 /**
  * render() contains the code necessary to render a convex hull in OpenGL.
  * @param (const ConvexHull<TYPENAME>&) obj
  */
 template <typename TYPENAME>
 inline void render( const ConvexHull<TYPENAME>& ch )
 {
    glColor3f( 1.0f, 0.0f, 0.0f );
    
    size_t numberOfTriangles = ch.triangles.size();
    for ( size_t i = 0x0; i < numberOfTriangles; ++i )
    {
        const Triangle<TYPENAME>& tri = ch.triangles[i];
        Render_Shape::Triangle_Solid( tri.getEndPointA(), tri.getEndPointB(), tri.getEndPointC() );
    }
    
    glColor3f( 1.0f, 1.0f, 1.0f );
    
    for ( size_t i = 0x0; i < numberOfTriangles; ++i )
    {
        const Triangle<TYPENAME>& tri = ch.triangles[i];
        Render_Shape::Triangle_WireFrame( tri.getEndPointA(), tri.getEndPointB(), tri.getEndPointC() );
    }
 }
 
 
 #endif // CONVEXHULL_H